Atomic Polarizability by Semi-Local and Non-Local Density- Functional Approximations

The static dipole polarizability of atoms with up to Z 1⁄4 20 electrons has been determined by density functional theory using prototypical local, semi-local and non-local approximations for the exchange-correlation energy. These successive generations of approximations provide very similar results for d; suggesting that their success in predicting atomization energies of condensed systems is mainly due to a better description of the atomic core region, while the description of valence sp electrons is largely unchanged. Density functional (DF) theory together with simple local and semilocal approximations for the exchange correlation energy Exc provides the basis for most of the computational studies of ground state properties of condensed matter [1]. Semi-local approximations adding gradient corrections to the popular local spin density approximation (LSD) for Exc predict cohesive energies and ground state geometries in fair agreement with experiments, as extensively documented in the literature [2]. The description of atomization energies, in particular, has improved dramatically in going from LSD to gradient corrected schemes, and there are clear indications of a further systematic (although quantitatively not very big) improvement in going from semi-local to fully non local approximations [3]. Despite the obvious success of these recent DF schemes, it is nevertheless meaningful to investigate whether the progress in atomization energies results from a better description of the highest energy occupied states (those directly involved in bonding), of the lowest energy atomic states, or of both core and valence states. A quantitative analysis of ground state properties other than atomization energies and geometries is required to answer these questions. We discuss here the static dipole polarizability ð dÞ of atoms, whose role is particularly important since this property is very sensitive to the energy, shape and symmetry of the electron states close to the Fermi energy, and polarization of the electronic cloud is the first elementary step towards the formation of a chemical bond. The sensitivity of d on the electronic configuration is emphasized by the large variations of this quantity along the periodic table, as shown in Fig. 1. Moreover, the static polarizability is related to optical properties by well known sum rules, and, for finite systems, it provides a measure of the spatial extension of the electron distribution [4]. We focus our attention on two recent approximations for Exc : (1) the semi-local Perdew–Burke–Ernzerhof approximation (PBE, see Ref. [2]), that is extensively used for applications, and (2) the meta-GGA approximation of Ref. [3], that introduces an explicit dependence of Exc on the local kinetic energy density of the occupied Kohn–Sham electrons. Polarizability is one of the first properties discussed in basic physics courses, and has a clear and intuitive definition. It is therefore surprising that dipole polarizability is not accurately known even for atoms. Most experimental determinations of polarizability rely on optical spectroscopy data. This route is very convenient for closed shell atoms, but far more challenging for open shell systems, resulting in large uncertainties in the experimental estimate of d: Similar difficulties affect computational schemes: the polarizability of closed shell atoms has been determined many times and by accurate methods, while open shell atoms have been far less investigated. The majority of the computational determinations of atomic polarizability relied on quantum chemistry many-body methods, including Hartree–Fock and post Hartree–Fock methods like configuration interactions, coupled clusters, etc. Hartree– Fock underestimates the polarizability of rare gas atoms, and overestimates it for the other elements. Post Hartree– Fock methods converge fairly rapidly to the (in principle exact) solution in the case of closed shell light atoms, while their accuracy degrades rapidly for open shell atoms or for heavy elements. This same trend is observed also for the computation of atomic polarizabilities. Systematic determinations of polarizability for closed shell atoms using DF methods have been pioneered by Mahan [5], Stott and Zaremba [6], Zangwill and Soven [7]. Numerical results were obtained mainly for closed shell atoms and ions, showing that LSD overestimates the atomic polarizability of rare gases by 10%; and underestimates the polarizability of alkali-earth metals. To the best of our knowledge, no systematic determination of d has been reported for the other elements. A non local recipe (self-interaction corrections to LSD, SIC) was implemented in Ref. [8], providing a systematic improvement in the evaluation of d; but, at the same time, manifesting fundamental and computational problems related to the SIC approximation. The application of e-mail: [email protected] Physica Scripta. Vol. T109, 166–169, 2004 Physica Scripta T109 # Physica Scripta 2004 more recent exchange-correlation approximations to the investigation of atomic polarizability has been very limited. Two representative studies are reported in Ref. [9] and Ref. [10], both based on the exchange-correlation potential of Ref. [11]. The lack of explicit computations is probably partly due to the difficulty of implementing efficient computational schemes suitable for the most promising but also most complex approximations, sometimes relying on implicit exchange-correlation functionals, or on orbital dependent potentials. The most accurate determinations of atomic polarizability by density functional methods have been performed using the Sternheimer approach [12], representing Kohn– Sham (KS) orbitals and their derivatives with respect to the external field on a radial logarithmic mesh centered on the atomic nucleus [5–6]. The resulting equations and their computational implementation are simple and elegant for closed shell spherical atoms. However, they become cumbersome for open shell atoms, and the representation of orbitals as radial functions times a single spherical harmonic implies that the unperturbed density is spherically symmetric. The method can be extended to nonspherical atoms, but its simplicity and computational efficiency are, to a large extent, lost. An additional difficulty is represented by the fact that the Sternheimer approach requires the analytic determination of the second functional derivative of the exchangecorrelation energy. This task is trivial for LDA, but it becomes increasingly difficult by including corrections of increasing sophistication. This problem becomes even more relevant in applying the Sternheimer approach to higher order response functions. The difficulty is not limited to the algebraic determination of the functional derivatives, but often includes numerical problems due to singularities in the exchange-correlation potential and of its functional derivatives. To overcome these difficulties, we implemented a program for the numerical minimization of the Kohn– Sham functional not relying on the analytic computation of any functional derivative of the exchange-correlation energy with respect to the density. The static dipole polarizability is computed from the dipole moment developed upon applying an an external field of the appropriate symmetry to the atom. The method is easily generalized to deal with perturbation of different symmetry (quadrupole, octupole, etc.) or to perturbations coupled to the spin density instead than to the electron density. Moreover, the computation of non-linear response coefficients could be performed by repeated computations of the induced moment as a function of the applied field. The method relies on the direct minimization of the Kohn–Sham functional for KS orbitals represented on a logarithmic grid for the radial coordinate, and on a discrete mesh for the ð ; ’Þ angular variables. The selection of the angular grid is performed in order to optimize the integration over the angular coordinates, as discussed in Ref. [13]. The computations presented below have been performed for grids of 2001 points for the radial coordinate, and 32 directions in the 1⁄20 1⁄20 ’ 2 domain. The angular mesh, together with appropriate weights, allows the exact numerical integration of all spherical harmonics Ylm with up to l 1⁄4 9 (see Ref. [14]). The combination of radial and of an angular grid allows to describe non spherical density distributions, either due to the perturbation or to a broken symmetry ground state. The core of the minimization algorithm is an efficient determination of the KS functional for any given set of N occupied orbitals